clear all;


n = 50000;

%s = RandStream('mt19937ar','seed', 3000);
%RandStream.setDefaultStream(s);

iota = 0 + 0.5*pi*rand(n,1);
psi  = 0 + pi*rand(n,1);
Phi0 = -pi + 2*pi*rand(n,1);

c = cos(iota);

a1 = -(1+c.^2).*cos(Phi0).*sin(2*psi)...
       - 2*c.*cos(2*psi).*sin(Phi0);


a2 = 2*c.*cos(Phi0).*sin(2*psi) ...
     + (1+c.^2).*cos(2*psi).*sin(Phi0);
 
 

a3 = -(1+c.^2).*cos(Phi0).*cos(2*psi)...
       + 2*c.*sin(2*psi).*sin(Phi0);

 
a4 = -2*c.*cos(Phi0).*cos(2*psi)...
       + (1+c.^2).*sin(2*psi).*sin(Phi0);
   
theorMean = sqrt(35)/4;

sqrt(mean((a2-a3).^2)) / theorMean 
sqrt(mean((a1-a4).^2))  / theorMean 
  
